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We present an efficient impurity solver for the dynamical mean-field theory (DMFT). It is based on 
the separation of bath degrees of freedom into the low energy and the high energy parts. The former 
is solved exactly using exact diagonalization and the latter is treated approximately using Green’s 
function equation of motion decoupling approximation. The two parts are combined coherently 
under the standard basis operator formalism. The impurity solver is applied to the Anderson 
impurity model and, combined with DMFT, to the one-band Hubbard model. Qualitative agreement 
is found with other well established methods. Some promising features and possible improvements 
of the present solver are discussed. 


I. INTRODUCTION 

In the past two decades, the dynamical mean field the¬ 
ory (DMFT) has been established as a powerful method 
for studying strongly correlated electron systemsfii^ In 
DMFT, a quantum lattice model for interacting electrons 
is mapped onto an effective single impurity model with 
a self-consistently determined bath. To study the multi¬ 
band lattice Hamiltonian or systems with strong spatial 
fluctuations such as two dimensional models, one needs to 
map the original Hamiltonian onto effective multi-band 
impurity modelsj^^— or onto a cluster Hamiltonian with 
many impuritiesFor DMFT and its extensions, the 
core is the solution of the impurity or cluster model for 
Green’s functions, or equivalently, for the self-energies<^ 
The efficiency of solving the quantum impurity or cluster 
models is essential for the applicability of DMFT and its 
extensions. 

Due to its obvious importance, fast and accurate 
method for solving impurity and cluster models has been 
the goal of intensive research activities since the devel¬ 
opment of DMFT. In this area, various quantum Monte 
Carlo (QMC) methods, such as Hirsch-Fye algorith m^°d^ 
and strongi^/weak-couplingi^ continuous time QMC are 
among the most powerful and widely used methods. 
They are numerically exact and can treat multi-band im¬ 
purity models and the multi-impurity cluster modelsJ^ 
However, these QMC methods cannot reach zero temper¬ 
ature and the dynamical correlations on real frequency 
axis are produced through numerical analytical contin¬ 
uation, which introduces additional errorsJ^ Numerical 
renormalization group methodic can produce accurate 
low energy properties for single impurity model, but the 
computational cost increase exponentially with the num¬ 
ber of orbitals or impuritiesiii The full exact diagonal¬ 
ization (ED)^ and the Lanczos methodic are relatively 
fast and easy to implement. However, at present, the 
total number of sites of the impurity/cluster model is 
constrained to about 9 for ED and 15 for Lanczos;^ This 
limits the energy resolution in the resulting spectral func¬ 
tion. 


Besides these more established methods, there have 
been attempts to develop new impurity/cluster solvers 
that are focused on the efficiency for complex impu¬ 
rity models, allowing for approximations. Among them 
are the weak coupling^! or strong coupling^^ perturba¬ 
tion theory, fluctuation-exchange approximation^^ non¬ 
crossing approximation^ Gutzwiller approximation)^^ 
double time Green’s function (GE) equation of motion 
(EOM) method)^^— superperturbation method;^ ex¬ 
tended recursion in operator space;^ and the distribution 
ED method,— etc. Recently, two methods appear to be 
able to go beyond the limitations in QMC, NRG and or¬ 
dinary ED. One is the ED algorithm based on natural 
orbitals and the other is the hierarchical equation of 
motion method— These two methods seem promising in 
DMET applications, but their implementations are com¬ 
plicated and the computational costs are relatively high. 

It is the purpose of this paper to present a new im¬ 
purity solver which is approximate but can be improved 
systematically. In this new method, the bath degrees 
of freedom are split into the low energy and the high 
energy parts. The influence of the low energy bath to 
the impurity is treated exactly using the traditional ED 
method. While the influence of the high energy bath is 
treated afterwards using an approximate GE EOM trun¬ 
cation scheme, that is, the alloy analogy approximation 
(AAA) — This idea is realized by using the standard ba¬ 
sis operator (SBO) and its GE EOM method^^ which 
is originally designed for extending the random phase 
approximation in the GF EOM formalism to the spin- 
1 Heisenberg model. Here we use it to derive a non- 
perturbative approximation on top of ED of a small sys¬ 
tem containing the impurity and Ug bath sites. We hope 
that this method can combine the advantage of ED, that 
is, high speed and being systematic, with that of GE 
EOM approach, that is, being flexible and giving con¬ 
tinuous spectral function. Our benchmark calculation 
shows that this method is fast and can produce qual¬ 
itatively correct spectral function using one bath site 
Us = 1. Therefore, this method is suitable for DMET 
study of more complicated models of strongly correlated 
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electrons, expecting reasonable results in a limited calcu¬ 
lation time. 

In the rest part of this paper, we first present the 
formalism in Section 11. The results of benchmarking 
calculations are presented and compared with ED and 
NRG results in Section III. In Section IV, we discuss the 
method and its possible extensions and applications. In 
Section V, a brief summary is given. 

II. FORMALISM 

In this part, we derive the formalism of the standard 
basis operator GF EOM method. For simplicity, we con¬ 
sider the single band Hubbard model on a Bethe lattice 
in the infinite coordination limit. The DMFT equation 
for this system can be obtained analytically and has been 
solved by many other established methods. It is an ideal 
system for benchmarking a new solver. 

The Hamiltonian of the single band Hubbard model on 
Bethe lattice reads 

H = d\^djc + Unj^mx - ( 1 ) 

i i,<y 

Here, d]^ and dia are the creation and the annihilation 
operators of a spin-cr electron on site i, respectively. The 
summation the nearest neighbour sites on 

a Bethe lattice with coordination number z. U is the 
on-site repulsion energy and /r is the chemical potential. 
In DMFT, this lattice Hamiltonian is mapped onto an 
effective Anderson impurity model of the form 

ddAim — ^ ^ T ^ ^ d^Cka^ 

k(7 k(7 

+Un^n^^ - E"- (2) 

<T 

Here Cka is the annihilation operator of the bath electron. 
The hybridization function is defined as 

A{u;) = J2VkS{<^ - ( 3 ) 

k 

with its counter part on the Matsubara frequency axis, 

r(*c.„) = ^—5^. (4) 

^ lUln - Cfc 
k 

A(aj) and are connected by the analytical contin¬ 

uation and the Hilbert transformation as below, 

A(u;) =- ImT —?> w -I- irj ), 

r(*a;„) = r ^^de. (5) 

J-oo - e 

We split the continuous bath degrees of freedom into two 
parts, a discrete part with 1 < k < Us and a continuous 


part with k > Us + 1. Here Ug is the number of discrete 
bath sites which will be treated by ED. Their respective 
hybridization functions are 

rig 

Ai(a;) = - €k), 

k=l 

A2(w) = A(w)-Ai(w). (6) 

Gorrespondingly, ri(ia;„) = Yl'k=i'^k / - (^k) and 

r 2 (iw„) = T{iLUn) - ri(iw„). 

Given a full continuous hybridization function A(a;), 
we first identify Ug discrete bath modes to define Ai(w) 
according to some importance criteria. The impurity 
problem composed of the impurity site and these Ug 
bath sites are diagonalized first. The influence of the 
residual continuous bath A 2 (a;) is regarded as a modifi¬ 
cation to the ED result and is treated by approximate 
method. This idea has been proposed by Hafermann 
et al.^ In their work, {I4,efc}(fc = l,2,...,ns) are ob¬ 
tained by a weighted fitting of the original hybridization 
function r(zu;„) by the rig-bath-site hybridization func¬ 
tion ri(ia;„). The influence of the residual hybridization 
T 2 {iuJn) is treated perturbatively, using skeleton diagram 
expansion up to third order in the interacting vertex in 
dual space. Although this approach produces rather ac¬ 
curate GF’s compared to QMG, the perturbative treat¬ 
ment of r 2 (*a;„) needs the two-particle GF X 1234 obtained 
from ED calculation as an input. For multi-band or 
multi-site impurity problems, the calculation and stor¬ 
age of X 1234 is expensive. This could hinder the practi¬ 
cal applications of this approach. In this work, we use 
GF EOM with truncation approximation, which involves 
only 2-time GF’s and thus avoids the direct calculation 
of 4-time GFs. For the parameterization of the discrete 
bath modes, we fit the low energy part of the full hy¬ 
bridization, which is important for the accurate descrip¬ 
tion of the Kondo resonance in the metallic phase (see 
below). 

In accordance with the splitting of bath, HAim is also 
split into two parts, 

dlAim — Hq Hi , 

ris ris 

Ho= + Z! Z! 

k—1 <7 k—1 <y 

+Un-fni - /ry^np., 

a 

Hi = 14 (cl^^da + dl^Cktr^ ■ 

fc>ns + l,(7 fc>ns + l,(T 

( 7 ) 

To integrate ED and GF EOM in a coherent way, we 
use the SBO formalism proposed by Haleyi^i^ Suppose 
Ho has been diagonalized and let and |/r) be the 
eigen energy and eigen state of Hq, respectively. That is, 
Ho\^) = SBO is defined as an projection operator 

in the Hilbert space of Hq in the |/i) basis, = \n){i'\. 
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Suppose that the Hilbert space of Hq is D dimensional, 
the total number of SBO’s is . The algebraic relations 
among these are 


The average ({H, 5}) can be determined self-consistently 
from the GF’s through the fluctuation-dissipation theo¬ 
rem, 


— A^fi 

A^i) Aoij^ = 5i)aAf^^ 

= ( 8 ) 

p- 

Clearly, the SBO’s defined above form a complete and 
linear-independent basis set in the operator space built 
on HqB Hilbert space. Any operator in this space can be 
expanded by Especially we have 

p 

( 9 ) 


Here the coefficients = {^jL\d„\i'). Usually the coef¬ 
ficients {/^y} form a sparse matrix and the sparseness 
can be used to simplify calculations. The superscript a 
in is used to remind us that this SBO appears in 
the expansion of da -. For the particle- and spin- conserv¬ 
ing Hamiltonian Hq, the excitation operators A^i, can 
be classified by its quantum numbers SN and SSz de¬ 


fined as 


N^A 


fll' 


= 6 NAfj_^ and 




fll' 


= 6 SzA 


flV ■ 


Here N and Sz are the total electron number operator 
and the z-component of total spin operator, respectively. 
The in Eq.Q are SBO’s with quantum numbers 

(diV = —1, 5Sz = —a 12). Here a = 1 denotes spin up and 
cr = — 1 spin down. In the following, we will use as 
defined here, and A^^ for general SBO’s with unspecified 
quantum numbers. 

In terms of SBO’s defined above, HAim can be written 


as 


EAim — Y^E.A E E 

fiy k—Ug+l 

oo 

+ E + dlcka^ , (10) 

k—ris + l 

where da = 'Zap and 4 = Writ¬ 

ten in this form, El Aim can be regarded as a new problem, 
describing a complicated impurity coupled to a reduced 
bath with the effective hybridization function T 2 (iu}n). 

The impurity Green’s function in Zubarev’s symbol is 
{{da\d)^))ui. It is expressed by the GF of SBO’s as 


{{da\dl))^ = EEW^ (11) 

7(5 

The EOM of a fermionic type double time GF ((A|H))^ 
reads 


( 12 ) 


1 /.oo / 

{BA) = -/ lm{{A\B))^+,r,du}. (13) 

^ J — oo ^ “r C 

Here 77 = 0-1- is an infinitesimal positive number. 

In the following, we apply the EOM to the SBO GF 
{{A'^lj\A^\))^. Appropriate truncation approximations 
will be introduced at certain levels to truncate the hier¬ 
archical EOM and produce a closed self-consistent equa¬ 
tions. The first order EOM reads 


( 14 ) 

To obtain the GF’s, we use the commutation relation 
[Aap^Ho) = {Ep — Ea) Aa 0 and the anti-commutation 
relation {Aap, As^} = dpsAaj + 5a-yAsp. Also, for 
fermionic type operators in different spaces, the follow¬ 
ing anti-commutation relations hold, Cka'} = 0 and 

^ka'} = 9. Using these relations, we get 


(uj -f Ea — Ep) ((A|)(^|A^J))^ — Si3s{Aa-y) + 5a^{Asp) 

fc>n+l,(7' 


(15) 


Here, {Aa-y) and {Asp) are thermal averages of SBO’s 
which conserve electron number and spin. The two new 
operators B^Z and D'^Z are defined as 


= {A-^, 4'} = E + E 

^ y 

Dap = {^a/3i da' } = E fpa^pP + E 

fi y 

(16) 

It is noted that B^p conserves N, while annihilates 
two electrons. Both are Grassmann even operators. 

On the right hand side of Ea. (fT5)) appear two higher 
order GF’s. {{B^p Cka'\A'^l))uj is the propagator of 
an electron from the small system Hq into the infinite 
residual bath, and {{D'^p cl.^,\A'^l))iz describes electron 
propagating inside the small system but being accom¬ 
panied by fluctuations of the residual bath. Both can 
be represented by GF’s of the form {{Afj_,yXka'\A'^l))ui, 
with A = c or A = cl. If one calculate EOM for 
these new GF’s, even higher order GF’s of the form 
((A^^Afco-'Apo.//|A^J)),^ will be produced. 

Different truncation schemes have been proposed for 
the Anderson impurity models under GF EOM formal¬ 
ism. AAA is made for GF’s involving one bath indexi^ 
The famous Lacroix approximation truncates the GF’s 
involving three bath index and gives qualitative descrip¬ 
tion for the Hondo resonancei^ Luo et al. proposed the 


u:{{A\B))^ = {{A,B}) + {{[A,H]\B)y 
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truncation scheme based on the concept of connected 
GF.— In their work the truncation at three bath index 
level produces improved results over Lacroix. 

In Ea. (fT51) . if we neglect the influence of the residual 
bath, i.e., if we remove the last two terms, the SBO GF’s 
are solved as 


{{AU\A^;s))u.= 


crt\\ _ ^/35(^a7) “b 


UJ ' 


■ En — Eft 


(17) 


Employing the fluctuation-dissipation theorem Ea. (fT51) 
and the SBO algebra Eq.®, one gets the self-consistent 
equations for the averages, {Aap) = /Zq. Zq = 

is the partition function of E[q. Putting the 
averages into the expression of SBO GF’s and using 
Ea. dTTI) . one gets 


term on the right hand side of Ea. dlQIl comes from the 
dynamics of the small system TLq accompanying the 
propagation of electrons. The EOM of this GF can 
lead to a chain of GF’s of the type {{OiCka'\A^\))uj 
(i = 1,2,...). This chain will get closed when the 
electron and spin conserving operator Oi traverse the 
whole space of Aap. The more advanced solution taking 
these into account will be studied later. Here, we make 
the decoupling approximation based on the argument 
that T 2 (iu}n) is supposed to be small and can be treated 
as a perturbation. 


(( 


DCr<7 TJ 

DCTCT tt 
^al3 ^^0 




{{drj\dl))uj 


1 {e-PE. + 

Zq u} + E^ — El, 

flU ^ 


(18) 


TDO-a- 


Hn 




= 0 . 


( 20 ) 


It is the Lehmann representation of the impurity GF of 
Hq, as expected. 

To take into account the influence of residual bath 
and improve over ED results, one needs to consider the 
EOM of the last two terms in Ea. (fT5]l . The EOM for 
^ka'\A'^l))ui , taking a symmetric form to keep the 
possible particle-hole symmetry, is obtained as. 




= (Kl,i?S^'c,.,}) + (( 


TDcrcr 


Hn 




-b 

-b 

-b 


Efc((TO', 












Ba'a ,da" 


Cfecr' j Cp^ 


(19) 


per" 


Since the most important effect of the impurity-bath 
interaction is supposed to be contained in Hq, we may 
devise simple truncation approximations for the residual 
bath. From the practical requirement of a fast impu¬ 
rity solver, we also need to avoid the complications in 
the EOM part. In the following, we introduce AAA ap¬ 
proximation which simplifies the EOM for spin-up GF 
by ignoring the commutator between the spin-down op¬ 
erators and the hybridization operators in H. That is, 
when studying the propagation of spin-tr electrons, the 
fluctuations in spin-d electron are frozen and treated as 
non-dynamical quantity. 

The last two terms in Ea. (fT^ come from 


(( 


DCTCr TT 


DtTtT TJ 




and ((Cfccr 

which describe the fluctuation effect from the impurity- 
residue bath hybridization when electron is propagating. 
In the spirit of AAA, they are ignored. The second 


For the first term in Eq. (fTUll ( J, Cka') , it could 

be calculated self-consistently from the GF of the type 
{{cka'\A'^^l))^, which is connected to ((A^^|A^J))„ using 
EOM. For the moment, to make the calculation sim¬ 
pler, we consider the following decoupling approxima¬ 
tion to this term, {A'^Ib^'^'C ka') « {A'^lcka'){B^'^') and 
{B'^P Cka' A"^^]) w -{B^p){A'^lckcr')- This leads to 


( 21 ) 

Finally, we ignore the spin flip process in Ea. ()19D since 
they are much less important than the spin-conserving 
propogation, {{B^p Cka'\A'^l))u^ ^ S^a'{{B^pCka\A'A^))ij. 
With all these approximations, the second order EOM is 
simplified into 

( 22 ) 

with = B^p. Carrying out similar EOM calcula¬ 
tion for the other GF in Eq.(IT5l), cl^,\A'^l))u, and 

making approximations alike, we obtain 

(23) 

with 

Putting Eg. d^^ and (1^51) into Ea. (fT51) . we obtain the 
closed set of equations for {{A'^p\A^\))ij as 

(w -b Ed — Ep) {{Adp\A^^g))= 6ps{Aaj) 

+ Saj{Asp) p ^ 

(24) 




























5 


Here, and are defined by = 

and {D‘^^ p,dl} = re¬ 

spectively. Their expressions are 




N'^a 

a.p,fii' 


I ra=¥ ra , ra* ra 

J uj3 J pot ' J otpJ 01 / 

(^f::rp}j +< 

I r(T* £(7 I £G* £G 

'Jl/0j pOL ' JotpJ 0U’ 


\ ^cr* £<7 

/ J J otT J pr 


\ ^ ra* £G 

/ ^ J rp J Tot 


(25) 


The coefiicients and satisfy two sets of 

exact relations which can be used for numerical test. For 
details see Appendix C. 

Multiplying to Ea. dMll and summing over 7 and 
S, we get 


{uj + Ea - Ep) {{A'^pld'l))^ — ^ f^p{Aa^) 

1 

+ T.f^s{Asp) + 

S py 

_ ^ (26) 


Eq. (|M)l is a set of linear equations about unknowns 
{((^a/sMo-))^}- For each w, it can be solved using it¬ 
erative methods. Parallel computation can be used for 
different w to reduce the computation time. It is the core 
equation on which our following calculations and discus¬ 
sions are based. 

There are unknown averages {{Aap)} in the inho¬ 
mogeneous term of Ea. (l26ll . They need to be solved self- 
consistently using GF’s via the fluctuation-dissipation 
theorem Ea. (fT^ . An equation about the averages can 
be obtained from each GF’s shown above, 


{d\Aip) = Y^rj{Asp) 

s 

1 r°° 1 

=-/ (27) 

^ J— 00 ^ 

Thus we obtain equations for the same number of un¬ 
known averages, which can be solved iteratively. In prin¬ 
ciple, from Eg. dMl) one could construct more than suffi¬ 
cient self-consistent equations for the averages {{Aap)}. 
However, it is noted that neither linear independence nor 
consistency is guaranteed in these linear equations, due 
to the approximations made for GF. There have been 
proposals for systematic ways of constructing sufficient 
and consistent equations;^ but in general this is still an 
open problem in the GF EOM approach. 

In order to keep our calculation scheme as simple as 
possible, we adopt the simplest approximation for calcu¬ 
lating the averages, 

{Ao^p) - {Ao^pf = (28) 


The advantage of using i7o result for {Aap) is that no 
self-consistent calculation is required and the sum rule is 
fulfilled exactly. The shortcoming is that the finite size 
effect of Ho will enter the final results and discontinuous 
change of averages and GF’s will occur when parame¬ 
ters such as ^ or C/ are tuned continuously at T = 0. 
For our calculation below, we stay at half-filling and this 
shortcoming has no effect. 

Putting the approximate averages into Ea. d^ . we get 
the linear equations for the final GF’s as 

= r.p [(^C.a)° + {App)°] , 

py 

( 29 ) 

where 

(^) — 

(uJ + Ect ^ 0 ) ^otp^0y 2 ^2 H" 2{^ ^)-^q 

Finally, the local GF of the impurity site is obtained as 

{{dM))^ = Y,r^p{{Aip\di))^. (31) 

OL0 

Before we present the numerical results, some discus¬ 
sions about the limiting cases are in order. First, in the 
limit Us = 0, equation Eg. (1^51) is equivalent to AAA, as 
being proved in Appendix A. This is because at Ug = 0, 
F 2 (w) = r(ci;) and the full impurity-bath hybridization 
is treated at AAA level. In the other limit Ug = 00 , the 
full continuous bath is contained in Hq and the residual 
hybridization r 2 (a;) = 0. As shown in Ea. dT^ . the exact 
GF obtains. We therefore expect a smooth interpolation 
between AAA and ED as Ug increases from zero. Espe¬ 
cially, the coherent quasi-particle features missing in the 
AAA spectral function should be produced for Ug > 1. 

Second, the set of equations Eo. d^ becomes exact in 
the limit U = 0, as proved in Appendix B. In the other 
limit of weak hybridization F(a;) = 0, our result is exact 
because the local impurity part is contained in Hq. 


III. RESULTS 


A. Anderson impurity model 


We first present the results of our impurity solver for 
the Anderson impurity model Eq.(2). For this model, we 
use a Lorentzian hybridization function 


A(a;) 


+ wg ’ 


(32) 


Wc = 1 sets the unit of energy and A describes the 
strength of the hybridization. The corresponding Mat- 
subara hybridization function reads 


T{i0Jn) 


TtAWc 

iuJn + iuJcSgn{ujn)' 


(33) 


/3./XJ/- 

(30) 
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In order to parametrize the discrete bath degrees of 
freedom in Hq^ we define the distance d between the in¬ 
put hybridization function r(iaj„) and that of the small 
system ri(zu;„), 

^ E - ri(*w„)| V<] , (34) 

n 

with N the number of summed Matsubara frequencies. 
d is then minimized with respect to variables and I 4 
{k = 1, 2,rig) use conjugate gradient method. For the 
number of bath sites Ug < 2, our results are sensitive to 
s value. We checked that s > 3 always leads to too small 
Cfe’s and below we use s = 2i^ 

Fig.l depicts the local density of states at T = 0.27rA 
and for different U’s, obtained using rig = 1 at the 
particle-hole symmetric point ^ At [/ = 0, the 

density of states coincides with the exact one, with the 
height of the peak p(0) = l/(A7r^). As U increases, the 
weight of the central peak transfers to higher energies, 
leading to a sharp central Kondo peak. The two broad 
peaks around U = ±U /2 are the incoherent Hubbard 
peaks, forming the well known three-peak structure. All 
the spectral function curves are non-negative and the sum 
rule J p{uj)duj = 1 is fulfilled. This shows that our ap¬ 
proach has no causality problem as in the strong coupling 
expansion.!^ 

As U increases, p{0) stay very close to the U = 0 
value. This is in qualitative agreement with the Fermi 
liquid behavior. At the quantitative level, however, the 
results in Fig.l deviate from NRG results (not shown 
here) in that p(0) decrease too slowly with increasing U 
at T = 0.27rA. This is because the Kondo resonance 
obtained here is mainly from the hybridization between 
the impurity and the single bath level contained in Hq. 
The whole low energy bath is only represented by this 
bath level at e = 0. The approximate way of calculating 
{Aaa) in Eg. d^ also makes our result less accurate in 
the temperature dependence. 

The quasi-particle weight of the impurity electron is 
defined as 


cIReECT(a;) 

-1 

Ir7iSc,(iu;o)’ 

[ duj \ 


Wo 


I]cr(iw„) is the Matsubara self-energy 
E.(*cu„) = [{{dM)lJ~' -[{{d.\di)h.^r\ (36) 

where {{da\d'l))igj^ is the impurity Matsubara Green’s 
function at t/ = 0. In Eg. (1551) . z is obtained either by 
evaluating the partial differentiation at the smallest pole 
of GF (ReflUld^), or from Matsubara self-energy E(iwo). 
Here we find that due to the finite energy resolution of 
small rig, the latter approach gives better result, as shown 
in Fig.2. z decreases monotonically from 1 at C/ = 0 
to zero in large U limit, in qualitative agreement with 
the NRG result. However, due to the same reason of 
small rig, our z value sensitively depends on T. Only at 



(0/(7tA) 

FIG. 1: The spectral function of Anderson impurity model at 
different U’s and temperature T — 0.27rA. They are obtained 
using ris = 1 for Lorentzian hybridization function with ttA = 
0.02. The broadening parameter is 77 = 10”"^. 


intermediate temperature T = O.IttA is our z close the 
NRG data which is obtained at T = 10“®. At low T our 
result is small than NRG result. 

We compare the spectral function p{uj) and the Mat¬ 
subara GF G{iujn) from different number of bath sites 
rig in Hq. The results are shown in Fig.3 at U = 37 rA. 
In Fig.3(a), we see that for rig = 0, there is no Kondo 
resonance at all because the full bath is treated at AAA 
level. For Ug = I, a sharp resonance appears at a; = 0. 
For rig = 2, the sharp resonance is replaced with a broad 
peak. NRG result at the same temperature is shown for 
comparison. Qualitative agreement between NRG and 
SBO-EOM results are observed, but systematic conver¬ 
gence with increasing rig is not obvious up to rig = 2 here. 
In Fig.3(b), Matsubara GF is compared with ED results 
using rig = 6 which is convergent already. It is seen 
that the SBO-EOM result achieves quantitative agree¬ 
ment with ED at rig = 1 level already. 

We also compared p{uj) with NRG data for larger U 
values. There, rig = 1 gives qualitatively good result, 
but for rig = 2, instead of Kondo resonance at a; = 0, 
two sharp peaks appear around some uip 7 ^ 0 , which is 
actually related to the position of poles of ri(w), i.e., e^’s 
in Hq from fitting. The poles in Fi(u;) will appear in the 
final results, by entering in Eo. (l30l) through 

r 2 (u;). In this sense, the sharp peak at w = 0 observed 
in Fig.l should also be closely related to the fact that 
Efe = 0 in Hq. Although the convergence of both p{ijj) and 
G{iuJn) in the large rig limit is guaranteed by the formula, 
qualitative correctness of p{uj) is limited to rig = 1 for the 
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FIG. 2: Quasi-particle weight z for the Anderson impurity 
model with Lorentzian hybridization function at ttA = 0.02. 
It is obtained using Us = 1 at T = O.IttA (sqnares) and 
T = O.OIttA (dots), respectively. The diamonds with guiding 
line is the NRG result at T = 10“®, obtained using K — 2 and 
keeping 256 states. 


moment. For rig = I, our calculation takes seconds on 
a PC. As Us increases, the computational cost increases 
exponentially and the efficiency advantage of the present 
solver will be lost very quickly. Therefore, we confine 
ourself to rig = 1 when we use this solver for DMFT 
applications. 


B. Hubbard model 

In this section, we combine the impurity solver with 
DMFT calculation to study the paramagnetic phase of 
the single band Hubbard model on the Bethe lattice. Due 
to the simplicity of the DMFT equation in this case, it 
is an ideal system for benchmarking our method. The 
Hamiltonian of Hubbard model reads, 

H = - tijcl^Cja + Un^ni - ^ ^ (37) 

(ij) ^ 

For Bethe lattice with infinite coordination number, the 
bare density of states reads, 

DW) = (38) 



CO 

n 


FIG. 3: The local density of states (a) and Matsubara Green’s 
function (b) of Anderson impurity model, for bath numbers 
Us — 0 (square), 1 (circle), and 2 (triangle), s = 2 is used 
in (b). Both are for Lorentzian hybridization at parameters 
[/ = SttA, T = 0.27rA, and ttA = 0.02. The solid line in 
(a) and the stars in (b) are NRG and ED results at same 
parameters, respectively. NRG result is obtained using full 
density matrix algorithm at A = 2 and keeping 256 states. 

Putting this into the standard DMFT self-consistent 
equations, it leads to the equation 

r(a;) = —G.(^). (39) 

The half bandwidth W is set as the energy unit W = 1. 

First we investigate the local density of states using 
rig = 1. Its evolution with U is shown in Fig.4. At 
[7 = 0, the exact p{uj) is obtained, as expected from Ap¬ 
pendix B. When we increase U, the weight of the local 
density of states at small frequency transfers to higher en¬ 
ergies, leading to a sharp middle peak, together with pro¬ 
nounced upper and lower Hubbard peaks at w ±C//2. 
This trend is consistent with previous studies using many 
well established methods, such as the iterative perturba¬ 
tion theory (IPT) and the quantum Monte Carlo (QMC) 
method)^ At the same time, p{uj = 0) declines even for 
U much smaller than the Mott transition point Uc- This 
is contradict to the Fermi liquid theory at such a low 
temperature. When U is close to 1.2, a gap opens with a 
small resonance peak standing inside the gap. The gap is 




























FIG. 4: The density of state of Hubbard model on Bethe 
lattice, obtained using SBO-EOM with Us = 1 as impurity 
solver. From top to bottom at a; = 0, 17 = 0.0, 0.6, 0.9, 1.2, 1.6 
and 2.0, respectively. Here fi = U/2 and T = 0.02, W = 1.0. 


fully formed at about Uc « 1.3FF, completing the Mott 
metal-insulator transition. 

It is seen that our results qualitatively resemble the 
ones of IPT and QMC. However, there are important 
quantitative difference. The obtained critical Uc « 1.3W 
in our calculation is much smaller than ~ 2.94FF 

from DMFT(NRG)^ and Ki 3.OFF from the two-site 
DMFT— . Compared to the AAA results = W, 

our result is between DMFT(NRG) and AAA values. 
This shows that using = 1, our solver is an improve¬ 
ment over AAA but is not as accurate as the variational 
method that employs one bath site in an optimal way^. 

Fig.5 shows p(0) as a function of U at T = 0.02FF 
which is much lower than the critical temperature. p(0) 
decreases monotonically and drop to zero at Uc « 1.3FF. 
The Mott transition in infinite spatial dimensions is of 
first order at finite temperatures. There is a finite regime 
Uci < U < Uc 2 where both metal and insulator coexist 
and there is metastable structure in thermodynamics4^ 
In the inset of Fig.5, we show that indeed our impurity 
solver can describe such first order phase transition, al¬ 
though the obtained coexisting regime is much smaller 
than what was found in previous studies. As U decreases 
across the transition point, p(0) jumps from zero to a 
finite value at different t/c’s, giving out a coexistence 
regime of U. 



FIG. 5: p(0) as the function of U at T — 0.02. Inset: the 
hysteresis near the Mott transition point. 


IV. DISCUSSIONS 

We first discuss the convergence of our results with re¬ 
spect to increasing n^, the number of bath sites in Hq. 
Gomparing p{uj) among = 0 to 2, it is found that 
apparently it does not show the expected convergence. 
This is mainly because the finite size effect in Hq will 
influence the final result, by adding poles in the residual 
hybridization function r 2 (a;). Therefore, for the real fre¬ 
quency spectral function, our solver produces an effective 
broadening for the ED results. On the other hand, the 
Matsubara GF converges quickly to the exact curve as Ug 
increases. Note that by construction, our method will be 
exact in the large Ug limit. The way of convergence can 
be understood from Eg. (130)1 : at those frequency points 
where the fitting is perfect ri(a;i) = F(a;i), our result 
equals to ED result because F 2 (aJi) = 0. For a given 
rig, it is empirically observed that there are Ug different 
(jji’s where r 2 (a;i) = 0, meaning Ug different frequencies 
where G{uJi) is identical to ED values. With increasing 
rig, G(w) thus approaches the ED result consistently. 

Second, the temperature effect is partly taken into ac¬ 
count in out calculation (not shown here). We have 
observed that the sharp Kondo peak for the Anderson 
impurity model disappears as T increases, leading to a 
broad featureless maximum. This is consistent with the 
expected behavior when T goes above the Kondo tem¬ 
perature Tfc. However, the hight of the spectral function 
/o(0) is higher than NRG result in this process. This 
problem, as well as the temperature dependence of z(U) 
curve can be partly traced back to the approximations 
used to simplify the self-consistent calculation of the av¬ 
erages (Aap), Ea. (|^ . We expect that by replacing this 
approximation with the full self-consistent solution, the 
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temperature dependence will be improved. 

Third, we discuss possible improvement that we can 
make in the future. In the approximations of GF, the 
spin flip process ignored in the derivation can be taken 
into account without additional numerical effort. On the 
computation side, our calculations for Ug = 2 takes only 
seconds on a PC. Increasing Ug will lead to much slower 
computation and larger storage requirement. However, 
we expect that the odd number of n = 3 and 5 can be 
reached with ease, which may give better description of 
both the Kondo resonance and the Hubbard peaks. The 
check of the convergence of results with respect to odd Ug 
will be done later. Also, with full self-consistent calcula¬ 
tion of the averages {Aap), we expect better description 
of the temperature dependence will be achieved. 


Putting these into Ea. (l25l) . we get the only non-zero el¬ 
ement among as = 2, and that among 

^ 24 ,tJ.v Af24 24 = 2. Also One obtains = 

^24,/.!/ = 0 for all (/r^)’s. 

With these information, Ea. (l26l) becomes 




a7 

w -f Ect — Ep — r(a;) 


(A3) 


Using the fact that A44 -|- A22 = and An + A33 = 
i-dldu one obtains the final GF at ris = 0 as 


((dtMj)). 


1 - W) 

Lu + pL — F(a;) 


(»t) 

(jj + pL — U — r(a;) 


(A4) 


V. SUMMARY 


summarizing the spin-f and -\. case, we obtain 


In this work, we develop and benchmark a new impu¬ 
rity solver for DMFT, based on the EOM of double time 
GF of SBO’s. Applying this method to the Anderson 
impurity model, we obtain qualitatively correct results 
for the real frequency spectral function p(w), the quasi¬ 
particle weight z, and the Matsubara GF. We applied this 
method to the one band Hubbard model on Bethe lattice 
through DMFT and obtained qualitative description for 
the Mott transition. Directions of further improvement 
in the future is discussed. 


{{d^\dl))^ 


1 - {ng) 

uj + fj, — F(a;) 


-f 


__ 

u! + p. — U — F(a;) 


(A5) 


This is exactly the expression of AAA if the average {ng) 
is calculated self-consistently. If we use its Hq value (ng-)° 
as was done in Ea. dMl) . it is not the full self-consistent 
AAA in general. But at exact half-filling (rig) takes its 
universal value 1/2, the resulting G{u!) recovers AAA 
again. 
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Appendix A: Solution for Ug = 0 case 

In this appendix, we prove that at rig = 0, Ea. (l2^ - 
m is equivalent to AAA . For rig = 0, Hq = Un-fni — 
/i^^ricr. Hence Fi(u;) = 0 and F2(aj) = r(w). The 4 
eigen states of Hq are denoted as 

|l) = |t), |2) = U), 

|3) = |0), |4) = |U). (Al) 

The corresponding eigen energies are Ei = E 2 = —p,, 
U3 = 0, and E 4 = U — 2pL. We decompose dg into SBO 
and obtain d'^ = A31 -|- A24 and ^4, = A32 — A14. The 
nonzero coefficients are 

III = fl = 1 , 

/l4 = /32 = 1- 

(A2) 


Appendix B: Exact solution at U = 0 limit 

In this appendix, we prove that Ea. (IM)l gives exact GF 
in the limit [7 = 0. 

In the case [7 = 0, Hq can be diagonalized on the 
single-particle level and we suppose that after diagonal- 
ization, Hq reads 

Hq = 'y ( CgQ,\^jQ,gg. ( 1 ^ 1 ) 

scr 

Here the operator Ugg is the annihilation operator of an 
electron in the molecular orbital s and with spin a. The 
impurity annihilation operator dg can be expanded as 
= J2g oisadsa, with Y,g \c(ga\'^ = 1- The eigen energies 
{U^l’s are the sum of the occupied single particle levels. 
The SBO Aap and the excitation energy E^ — Ep are 
defined similarly as in Eq.Q. In this case, considering 
that [agg, Hq] = Cggagg, the expansion of Ugg in terms of 
SBO’s {Alp) reads 

= (B2) 

Here denotes the summation with the constraint 

Ev — E^ = eg, i.e., the summation of degenerate excita¬ 
tions only. 
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Multiplying Eq. (1261) with and summing over a and 
j3 under the constraint Ep — = Cs, one obtains 


(w - €s) {{asa\dl))u, 

= {{a-scr^dl}) + ^T2{u}){{{{asc,,di},da}\dl))^ 

-ir2(-w)(({{a..,d.},4}|4)).. (B3) 


E d.} = {{d., 4}, d4 = 2d., 

q:/3 

T.fp:{Kp,dU = {{4,d4,4} = 24, (Cl) 


q:/3 


Here we have used the relation jg easy to obtain the following relations 

{dd^p^drj} = dj.}, d.}, and “ 


{Kp^4} = {{Kp^da},di}. 

Using the definition of Usa, Ea.(IB3ll leads to the equa¬ 
tion for GF as 

E] fapd^ap,^iu — 24i^, 
al3 


{(rf<rl4))^ = (e!^) [l + r2(u;)((d.|4)E . 

a/3 

(C2) 


(B4) 

Note that for [/ = 0 case, the GF of the small system EIq Similarly, from the operator equations 
can be expressed as 


Go(a;) = E 


U! — €s U! + ^ — Ti{uj) 


(B5) 


Putting this equation into Ea. (IB4l) . we Hnally recover the 
exact impurity GF of El Mm at U = 0, 


{{da\dl))^ = — 


1 


1 


Go^(u;)-r2(w) uj + fi-T{uj) 


■ (B6) 


Appendix C: Some Exact Relations 

From their definitions, we can derive the following ex¬ 
act relation about and Nap^^v From the opera- 


E fp*c.{Kp. da} = {{ 4 , 4 }, da} = 0 , 

E 4 } = {{d., da}, 4 } = 0, (C3) 

one can get the relations 

E! fpa^d^oiP,fiv = 0, 

a/3 

E! faP^aP,fj.u = 0. (C4) 

a/3 
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